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Abstract 

The Hybrid Monte Carlo algorithm is adapted to the simulation of a system of classical degrees 
of freedom coupled to non self-interacting lattices fermions. The diagonalization of the Hamiltonian 
matrix is avoided by introducing a path- integral formulation of the problem, in d-f 1 Euclidean space- 
time. A perfect action formulation allows to work on the continuum euclidean time, without need for 
a Trotter-Suzuki extrapolation. To demonstrate the feasibility of the method we study the Double 
Exchange Model in three dimensions. The complexity of the algorithm grows only as the system 
volume, allowing to simulate in lattices as large as 16^ on a personal computer. We conclude that 
the second order paramagnetic-ferromagnetic phase transition of Double Exchange Materials close to 
half-filling belongs to the Universality Class of the three-dimensional classical Heisenberg model. 

05.10.Ln, 75.10.-b, 75.30.Et 

1 Introduction 

Most of the models so far proposed to study the Colossal Magnetoresistance manganites (CMR) Q, 
share an extremely simplifying feature: an assembly of non self-interacting lattice fermion is coupled to 
an extensive number of classical continuous degrees of freedom (the localized core spins of the Kondo 
model and of the Double Exchange Model and/or the Jahn- Teller lattice distortion fields j^). Other 
physical context where this simplifying feature appear are the pyrochlores or doubles perovskites. 

The non self-interacting nature of the electrons in these models makes it possible to explicitly perform 
the trace in Fock space, in terms of the single-particle eigenstates. This yields a positive Boltzmann 
vi^eight for the continuous classical degrees of freedom, that for the sake of brevity we will call spins 
in what follows (although they could be a lattice distortion field!). In principle, the resulting problem 
could be simulated by means of a Metropolis algorithm. However, the update of a single spin requires a 
diagonalization of the single-particle Hamiltonian matrix, which has a computational cost proportional 
to the square of the lattice volume (if the most sophisticated available algorithm is used). This implies 
that the time needed to update all the spins on the lattice scales at best with the cube of the lattice 
size. This problem has prevented the study of systems with more than (say) two-hundred spins (a 6^ 
lattice) in the simplest of the above quoted models, the Double Exchange model (DEM), although most 
simulations ^ are done with a hundred or less spins. This is certainly not enough for an accurate 
study of phase-transitions where most of the interesting physics occurs. 

In this paper, we reformulate the problem in the path-integral formalism, obtaining an exact represen- 
tation on d-t-l dimensions for the fermions and d dimensions for the (classical) spins. In this representation 
a positive Boltzmann weight is obtained, and the update of the spins can be done by means of the Hybrid 
Monte Carlo (HMC) algorithm |^. For the DEM, the computational cost of a full-lattice updating is 
empirically found to grow as the lattice volume (although a worst-case estimate would have yielded a 
square- volume growing). In addition, the autocorrelation time for HMC is proportional to the correlation 
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length while with the Metropolis algorithm in the Hamiltonian formalism it grows like the correlation 
length squared. We will show that a standard simulation on a 4"^ lattice yields fully compatible results 
with our HMC algorithm, but the latter allows to simulate a 16^ lattice on a personal computer. In this 
way we are able to obtain meaningful results for the phase diagram of the DEM model. Some attention 
will be paid to the largeness of the finite-size corrections on the small lattices. We will also show that 
in the absence of superexchange coupling between the spins (whose numerical treatment is straightfor- 
ward), the Double Exchange Model near half- filling presents a second order phase transition between 
the paramagnetic (PM) and the ferromagnetic (FM) phase, that belongs to the Universality Class of the 
three-dimensional Heisenberg model. Work is in progress for the study of the phase-diagram of the DEM 
complemented with a first-neighbors antiferromagnetic superexchange interaction. The final goal is to 
confirm that the antiferromagnetic coupling is able to turn this PM-FM phase transition from second to 
first order as predicted by Mean-Field |^ . The phenomenological importance of reliably finding PM-FM 
first-order transitions between phases of very similar electronic densities cannot be overemphasized . 

The structure of the paper is as follows. In the next section we describe the DEM, introducing our 
notational conventions and deriving it from the Kondo lattice model. This somehow academic exercise 
will allow to introduce in a natural way a mathematically equivalent formulation of the DEM in terms of 
SU(2) matrices rather than in terms of classical fixed-length spins (to this respect, appendix^ will be also 
of interest). This representation of the model will allow for an enormous improvement of the numerical 
stability of the integration of the equations of motion during the Molecular Dynamics part of the HMC 
algorithm. In section 3, we present the path-integral formulation of the model, and prove its mathematical 
equivalence with the Hamiltonian one. In section 4, we give details of our implementation of HMC. Section 
5 is devoted to consistency checks: we show numerically how our perfect action formulation avoids the 
need for a Trotter-Suzuki extrapolation to continuum Euclidean-time and we compare the numerical 
results of the HMC simulation with an usual Hamiltonian one. In section 6 we present our results for 
the PM-FM phase transition at half filling. Section 7 is devoted to conclusions. We also include three 
appendices with useful formulae and the proofs of some relations used. 



2 The model 

We consider the lattice Kondo model on a cubic lattice of side L and volume V ~ L^, where periodic 
boundary conditions are applied. On each lattice site we have a classical localized spin, (pg. of unit-length. 
The spins interact with a band of lattice fermions through the Hamiltonian 

^ = ^^<^i,aHx..a;y,l3Cy,f3 , (1) 
X.a y.p 

where x and y run over all nodes of the spatial lattice, and a, /3 = 1, 2 are spin indices. The single-particle 
Hamiltonian matrix consist of a hopping term plus the Hund coupling with the localized spins: 

d 

Hx,a-y,p = —t'^^5a,l3 [5x:.y+l + 5x:y-i\ — Jh Sx;y{(f'x ' (^)a,f3 , (2) 
i=l 

where i is the unit vector in the i direction and a = {cti, 12, cr^} are the Pauli matrices. For the particular 
case of the CMR manganites, the localized spins represent the three core manganese t2g electrons that, 
due to the Hund rule, yield a = 3/2 spin that for most purposes can be considered as classical. The 
conduction electrons, represented by the creation and annihilation operators Cx,cn cj, ^, occupy the lowest 
of the two manganese eg orbitals, split by a Jahn- Teller distortion. 

The statistical properties of the system with an explicit superexchange antiferromagnetic coupling 
between the localized spins can be obtained through the partition function. Choosing units such that 
fee = li the partition function reads 

'D^e-T^'^TrFockg-^fw-MA^)^ 
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the superexchange hamiltonian being 



''X ^x-\-t 1 



(4) 



and M is the number operator 



W,n] = 0. 



(5) 



The problem can be enormously simplified, due to the non self-interacting nature of the Hamiltonian 
(|l|). Although 7i is a 4^ X 4^ matrix in the Fock space, the trace in Eq. (^) can be explicitly taken if the 
eigenvalues of the 2V x 2V single-particle Hamiltonian matrix, {En}n=i....,2V i a-re known: 



Tr 



Fock g_i(W_^AA) 



i:„ log 1+' 



(6) 



It is thus clear that, as we have said in the introduction, the resulting Boltzmann-weight is positive, and 
that the model can readily be simulated by the Metropolis algorithm, up to the computational caveats 
mentioned in the previous section. For numerical calculations based on this strategy, see Ref. 

The dimensionality of the matrices can be still reduced in a factor of two, in the limit of large Hund 
coupling, thus obtaining Zener's double-exchange model One first makes a unitary transformation 
that diagonalizes the Hund coupling term in Eq. (0): 



H 



x,a;yl3 — ^x.yU (^4>x)a,(i 

U{$) = 



COS -e'(-+^)/2 
2 



sin-e'(-+^)/2 
2 



sin-e'^'^-^)/^ 
2 



.cos-e'(^-^)/2 
2 



(7) 
(8) 

(9) 



where and ip are respectively the polar and azimuthal angle that determine the spin direction. It will 
also be important in what follows our choosing of U{(j)) as an SU(2) matrix. The resulting single-particle 
Hamiltonian matrix is 



r 

Hx,a;y,0 = - Jh {0-3)0,13 -t'Yl (U{(fx)U''{$y)) 5x.y+i + (U{(j)x)U\$y)) 5^,y^i 



(10) 



Due to the largeness of the Hund coupling one should keep only the electron state with spin parallel to 
its core spin (the "1" state in the representation of Eq. (|o|)). The truncated single-particle Hamiltonian 
matrix is then 



a 

Hx,y = ~tJ2 (c^(4)t/^(4-i))^/.,y+i + (c/(4)f/^(0 



x,y — % 



(11) 



Let us take a look at the product 
(C/(0,)f/t(0^))^^ 



cos ^ cos ^ + sin^sin^e-'(^--^«) 
2 2 2 2 



(12) 



The term between square brackets is nothing but the hopping term of the DEM model (see e.g. Q). 
Thus we see that the matrix in Eq. ( pT| ) is actually an unitary-transformed of the usual hopping term, 
the unitary transformation being 

(13) 



^xy ^x,y^^^^ 



Now, the expression in Eq. (y_l|) is extremely more convenient for an HMC study than the usual one. 
In fact, during the Molecular Dynamics part of the algorithm, one needs to take care of the constraint 
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{4>x)^ = 1. It can be done with a modification of the usual equation of motions as shown in Ref. To 
get these new equations of motion one needs to express the hopping term of the DEM in terms of the 
Cartesian coordinates of the spin ((/)^, (lP,<j)^) 

cos^cos^ + sin^sin^e-'("v)^;l I ^/l + cj^U/l + + '^"^ , ' \ ■ (14) 

Indeed, a working HMC algorithm can be obtained using the above representation which is not 
analytic at the sphere South Pole. However, during the Molecular Dynamics step of the HMC, one 




needs the derivatives of the right-hand side of Eq. (14), which at the South Pole are even more singular 



than (14), resulting on a poor numerical stability of the integration of the equation of motion. On the 
contrary, the expression of the hopping term as a function of the SU(2) matrices is smooth. Moreover, 
as discussed in Appendix nothing changes if we substitute the integrations over the spin-field in the 
partition function, by an integration over the SU(2) group. If needed, the spins 4)x can be obtained from 
the SU(2) matrices using the formula (see appendix ^) 

'/'i = ^Tr(a,;7tfT3C/.) , J = 1,2,3. (15) 

Thus we will consider the following statistical system, which is strictly equivalent to Eq. (^ in the 
double-exchange limit: 



r -iffS^+i:^ , log i+e — "T 

Z = DUc ^ V J , (16) 

^SE ^ :^^^Tr[(C/ta{/.).(L/I^,aC/.+,)] . (17) 

In the above expression, T is the temperature and En are the eigenvalues of the single-particle Hamiltonian 
matrix defined in Eq. ([Tl|). 

Although the SU(2) field Ux is still a constrained variable, it can be dealt with using well established 
techniques from lattice-gauge theory jiot . 

Let us also finally mention that the single-particle Hamiltonian Eq. (|l^), is unitary equivalent to 
minus itself, the unitary transformation being (x, y, z are the lattice coordinates of a;) 

Ux.y = 5x.,y{-ir+y+' (18) 

This ensures that the spectrum is symmetric around zero and therefore half-filling corresponds to /i = 0. 



3 From the Hamiltonian to the Path-Integral formulation 



In this section, we will show how to obtain a numerically tractable path integral representation of the 
partition function (|l^) although our results will be valid for the general problem outlined in the introduc- 
tion: classical continuous degrees of freedom coupled to non self-interacting fermions. In subsection [3.l| 
we shall also explain how some important fermionic observables can be recovered in this formalism. 

Let us first state the following well known expression for the partition function in terms of a pair of 
anticonmuting Grassman fields {\E'^(t), ^I^a, (r)}, where r is the Euclidean time |T3], 



Z = 

Sp = 



DU D-^ D-^^e-^^-T" 

T 

dr 



(19) 
(20) 



In the above expressions, H is the single-particle matrix defined in Eq. (|l l|) , and the Grassman fields 
verify antiperiodic boundary conditions in the Euclidean-time direction 



(21) 
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Now, in order to transform the representation (20) onto a numerically tractable problem, we introduce 



a time discretization. Wc introduce Lt time slices (for technical reasons, Lt will always be an even 
number, see Eq. (|3^)), with a spacing such that 

h 

Lrttr = — . (22) 

In this way, instead of a three-dimensional lattice and a continuum time, we have a four dimensional 
lattice. The Grassman fields now depend on the discrete coordinate Xr, 

T = XrCLr , X,- = 0, 1 , . . . , — 1 , (23) 

and verify the boundary conditions 

*L = -*L,. ' *o,. =-*L.,.. (24) 

The fields Ug. instead, being classical, do not depend on Euclidean time. In order to check how close our 
time discretization is from the continuous limit of Eq. ( pO| ) , we need to compare a^- with the natural time 
unit of our problem, h/t (see Eq. (pi])). Therefore, the dimensionless parameter that controls how close 
we are to the continuum-time limit is 

A = ^. (25) 

Our discretization should be such that in the A ^ limit Eq. ( pO| ) is recovered, much in the spirit of 
the Trotter-Suzuki extrapolation. From now on let us also adopt the convention that the quantities with 
dimension of energy, T, yU, Jaf and the matrix H are measured in units of t, in such a way that for instance 

With all our notational conventions settled, the discretized form of the action is 

- E e^'^l.-^-^+i- - E ^-..=»[exp(Ai?)],;^^..,„ (27) 

= T.T.'^L=^^^Lx-y.,y'^y^^y (28) 

Xt,x y-r,y 

The last equality on the above expressions defines the so called fermionic matrix M^. The rationale for 
including the chemical potential on the temporal link that joins the (xr, x) site with the {xr + 1, x) one, 
can be found in Refs. |Q. It is easy to check that in the A ^ limit the continuum-time action is 
recovered. The exponential form in the spatial part of Eq. ( p7| ) is preferred over more straightforward 
ones because it yields a perfect action, as shown below, without any time discretization effect. For the 
particular case of the DEM model, the action in Eq. ( p7| ) can be directly simulated. For other models, 
the approximated form 

[eMm]^.^ ^ 5^.,y + \H^,y, (29) 

could be the only feasible one, but it would makes mandatory to consider the A — s- extrapolation. 

To show the correctness of our path-integral, it is useful to first introduce the time Fourier transformed 
field, 

= -^E^'^'^'^^f^^- ' (30) 
where the sum extends over the Matsubara frequencies (see Eq. (p^)), 

2^ Lr-l 1 1 Lr-l 

p' = rj' '^^ — 2-'---'"2'2"--'^-- 

The fermionic action now reads 

5^ = E ^''^"""KoAux E *Po,. i^M^H)]^.,y %o,y ■ (32) 
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Therefore, the fermionic matrix defined in Eq. ( pSj) is block-diagonal in this basis 

Mi,.;p[,,y = V,Pi^»=;i/(Po,A), (33) 
A^;y{po,\) - e'^^+^P0 6^.y - [eMmUy, (34) 
and the Hamiltonian matrix, being Hermitian, verifies 

^t(po,A)=A(-po,A). (35) 
The partition function is then (using the Grassman version of Gaussian integration) 

= /]JdC/^ ]J d*,o^j, ]^ Ll^t^ _,e"^"^''o,.2,o.„<,..A/po,„;™ (3g) 

X qo,y po,x 

= I DU det [M^] e-^"'^ ^ /^^^ H '^^^ [A\pq, X)A{po, X)] . (37) 

'' Po>0 

It is then clear that the block-diagonal from of the fermionic matrix yields a positive-definite Boltzmann 
weight. Now, in order to relate this Boltzmann weight and our target expressed in Eq. (|6|), let us first 
notice that the eigenvalues of the A{pq), in terms of the V eigenvalues of the Hamiltonian single-particle 
matrix (|Tl|), are 

E^{X,Po) = e^'^+'P° - e^^" . (38) 



Now using the equation 



n 



(e^^A+.po_e>^) ^e'°H'+°' ' ^ (39) 

proved in appendix ^ we find for the fermionic determinant 

det[M^] = JJJKe^^+'Po -e^^") , (40) 



e^"W^+'="^)+-l, (41) 



^4^Trff E „ [log l+c- j I 



(42) 

Since the single-particle Hamiltonian matrix (pi] ) is traceless, it is clear that the discretized action exactly 
reproduces the target Boltzmann-weig ht (|l^ and thus it can be rightly called a perfect action. In the 
general case, though, one would have to take out by hand the TrH/T from the Boltzmann weight. 

Before ending-up, let us say a few words about the (in principle) non local matrix exp[AiJ]. In a 
model such as the DEM, where the eigenvalues of the single-particle Hamiltonian matrix are within some 
a priori known bounds, it can be numerically computed with a polynomial expansion as described in 
appendix |^. In other cases, although the energy must be bounded from below if the system is to be 
stable, an upper bound may not be available. Then one will be enforced to use an approximation such as 
Eq. (|2^) . It would anyway be convenient to add a multiple of the identity matrix to the single-particle 
Hamiltonian, in order to have a positive spectrum. From the above analysis, it follows that with the 



approximation Eq. (29), the simulation would be exact for the A dependent single-particle Hamiltonian 



H^^\ log(l + XH) . (43) 
A 

Therefore, there would be a deformation of the spectrum (as one finds using the Trotter- Suzuki formula 
at finite time-slicing), that would disappear on the A — *■ limit. More bothersome, there would be an 
empty-band dynamical effect. Indeed, even in the /i —oo limit, where fermions should not influence 
the classical degrees of freedom, the TiH^ term of Eq. ( ^2| ) would be present, and as one has 

ilog(l + Ai/) = i/-^i/2 + ... , (44) 
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the spins would definitively feel this spurious interaction, even if Tr_ff = 0. One can completely cure this 
problem, by using the Boltzmann-wcight 



Z= / L>C/ Jl det 



Po>0 



^t(pQ^ A, ^ = -oo)v4(po, A, ^ = -oo) 



e — — ^ (45) 



that can be simulated, using a straightforward modification of the HMC algorithm explained in section ^, 
because the matrices ^(po,/j,) commute for all values of po, A and /i. 

3.1 Fermionic Operators 
3.1.1 Charge density 

Let us call |a;) the state localized on the lattice site a;, and |n) the eigenvector of the single-particle 
Hamiltonian matrix corresponding to the eigenvalue The charge density on site a;, for the given 
configuration of the spin-field is 



p,^Y.\(^n\x)f^—^ , (46) 

e^J^ + 1 



n=l 



while the average charge-density on the lattice is 



p^\llp^-\Y.^^r-.- (47) 

Now, using (see appendix ^) 



r, ^ ouA-l-ipo _ pAB ^ -, , (48) 



we obtain 

P^ = tX ^ eviif-eAi^„ = ]^ E [^(^0, A)]- , (49) 

while 

p= -i-l^e^^+'«'Tr[^(po,A)]-i . (50) 

Po 

In practice, we make use of the equality between the thermal average of p-^ and the one of p, since most 
of computer-time during a HMC simulation is spent in the inversion of the matrices ^(pojA), which is 
done row-by-row. We therefore only calculate one row of the inverse matrix, and store the corresponding 
value of Pa, . 

3.1.2 Fermionic Energy 

In the Hamiltonian formalism the energy (per spin) for a given configuration of the spin degrees of 
freedom, is obtained from the logarithmic derivative with respect to the inverse temperature, /3, of the 
partition function, and has the form 

e = -y f"~^ (51) 



Using again Eq. (48), we can write the first term of the RHS of the previous equation as 

_}_ gAM+ipo 

Y /-^ /-^ gpA+ipo _ gAK 



and thus 



1 ^e^'^+'Po lTr[(i7-M)A-i(po,A)] . (53) 



Po 
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As in the case of the density, one cannot afford to calculate the full trace, but rely on the translational 
invariance and calculate 

eAx) = -l^e^^+'^'° [iH-f,)A-\po,\)],^, , (54) 

Po 

that can be readily obtained once the x^^ row of the matrix A{pq, A) is known from the density calculation. 

On the other hand, the total specific heat cannot be calculated in a practical way from the thermal 
fluctuation of the energy. Indeed, one can easily find that 

A representation analogous to Eq. ( p^ can be readily obtained for de-p/dj3. The real problem is the 
calculation of (e^), because we do not know ep, but e-p{x). It is easy to convince oneself that to substitute 
ep by ep^x) on the calculation of (e^) produces a systematic overestimation, magnified by the V prefactor. 



4 Our implementation of HMC 



In this section, we will give the necessary details about our implementation of the HMC algorithm 
The reader interested in a full exposition of the algorithm may consult jl^ . 
Let us recall that we want to simulate the statistical system 



I DU \{ det[At(po,A)^(po,A)] e""^ 



(56) 



Po>0 



As usual, the first step is to get rid of the fermionic determinant by using Gaussian integration, introducing 
the Lt/2 pseudofermionic (commuting) fields, ^PpQ,x'- 



det Af 



PO,!C 

\P0>0,!D / \po>0,!C / 



exp ■ 



^Po,a=(^^(Po,A)A(po,A))^>p„j, i . (57) 

po>0,x,y ) 



In our case of constrained variables belonging to the SU(2) group, one introduces 3V momenta (one per 
group generator and per lattice site |p^ ), by multiplying Eq. ( ^7|) by unity written in the form 



rr dP^ 



p2 
^ X 

' 2 



(58) 



So we end-up with a classical-mechanics model, that can be studied using the Molecular Dynamics 
method. In this model the kinetic energy is 



'1 

X I 



(59) 



while the potential one is 



U 



H 



SE 



T 



(60) 



po>0,x,y 



Following the standard procedure, at the beginning of each Molecular Dynamics trajectory, the momenta 
are extracted with the corresponding Gaussian probability (|5^), while the pseudofermions are obtained 
from a Gaussian vector S,pa,x as 

<Ppo = ^^(Po, A)^po . (61) 
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In practice, the pseudofermions being instantaneously thermalized, they are not changed during the 
trajectory. It is useful to consider instead two molecular-dynamics time dependent fields 

= (AtA)-V, (62) 
^ ^ Aij, (63) 

Although if is not changing, the matrix A{pq,X) changes when the field Ux follows the dynamic. The 
equations of motion adapted to the SU(2) group constraints are (the dxj derivative is defined in 
appendix ^) . 

Ux = (iPx • <?) C/cc , (64) 
Px,j = -dx,,U. (65) 

The hard part to calculate is of course 

(66) 

= -^\dx,,A^)i - i{dx.jA)r^ - -2Re [^^(9,,^- A)ry] 

Thus we see that a knowledge of the full inverse A{pq^ A) matrices is useless, and it is enough to consider 
the field r/ defined in Eq. (|6^). Once we know how to calculate derivatives of the exponential of the 
single-particle matrix (see appendix^), the rest of the calculation is standard: we numerically integrate 
the equations of motion by means of the SU(2) leap-frog algorithm inverting the A(po, A) matrices 
using a conjugate-gradient method. A numerical trick of some relevance is that one can calculate the 
inverses during the Molecular-Dynamics steps of the algorithm with far less accuracy than during the 
Monte Carlo accept-reject step ^, For the exponential of the single-particle Hamiltonian, we have 
used an order of the polynomial expansion such that the error is smaller than 2 x 10^"' all over the 
spectrum. 

An important remark about the algorithm is that pseudofermionic (four dimensional) variables 
can be straightforwardly generated following the exact probability distribution. This allows to simulate 
systems very long in the time direction without compromising the autocorrelation time. 



5 Consistency checks 

We consider in this section some tests performed to check the algorithm. Firstly, we should mention that 
although the computer code for the HMC is rather complex, most of the routines are very easy to check. 
For instance, the matrix inversion is self consistent, and the integration of the equations of motion can 
be directly checked as they should conserve the Molecular Dynamics Hamiltonian (up to second order in 
the leap-frog step). In addition, we have checked explicitly the reversibility of the equations of motion. 
A posteriori, it is very useful to control the Crcutz parameter p7| defined as 

{e-^""°) , (67) 

where all variations of i?MD must be considered (accepted or not). This quantity should be 1, and its 
measure is a very strong check of the simulation. A deviation would mean a reversibility problem or a 
lack of equilibration. We have readily checked this parameter in all the simulations. 

Regarding the comparison of the time discretized model with the physical continuous limit target, 
we have performed the following two types of test. On the first place, we have simulated a 4^ lattice at 
T = 1/8, for decreasing values of A, using Eq.(^9|), with a shift of the identity 6A that ensures a positive 
spectrum. The empty-band dynamical effect is avoided using Eq.(^5|). We have chosen /i = —3.5 and 
Jaf = which, for T — 1/8 is near the Paramagnetic- Ferromagnetic transition. The results are displayed 
in Fig. |l| for several quantities. For this selection of the parameters a linear behavior in A is observed 
only for large values of Lr = (TA)~^. We have also carried out a simulation with the Perfect Action (see 
appendix ^ in a 4^^ x 16 lattice at A = 0.25 with a 6 degree polynomial approximation. The result is 
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Figure 1: Continuous time limit for the nearest neighbor correlation, magnetization squared and charge 
density in a 4^ lattice at Jaf/^ — 0, T — 1/8, fi = —0.35. The left-most open symbols require Lr = 2048. 



The filled symbols correspond to a simulation with the Perfect Action (109) using A — 0.25, Lr — 32 



JAp/t 




-0.01 


0.05 


0.2 


0.3 


{Sx • Sx+i) 


Hamiltonian 


0.7734(8) 


0.3817(18) 


-0.4699(5) 


-0.6842(5) 




HMC 


0.7717(10) 


0.3838(8) 


-0.4697(3) 


-0.6852(4) 


{Mf 


Hamiltonian 


0.7149(16) 


0.0162(7) 


0.0130(4) 


0.3580(13) 




HMC 


0.7127(16) 


0.0152(3) 


0.01340(16) 


0.3611(11) 



Table 1: Comparison of the results of Hamiltonian and HMC simulations in lattice with T ~ 1/20 
at half filling (/i — 0). We show the correlation between nearest neighbor spins and the square of the 
magnetization (magnetization staggered when that correlation is negative). The numbers correspond to 
10000 measures in each case. 

plotted as filled symbols in Fig. |l]. The agreement is excellent. The selection of A for a Perfect Action 
simulation should be taken looking at the performance of the algorithm. Most of the results presented 
in this article have been obtained with A = 0.125 and a polynomial degree of 6. Larger values of A have 
the advantage of requiring smaller values of but the matrix inversion is more expensive. Conversely, 
smaller values of A require larger Lt while the benefit in the matrix inversion is scarce. 

Our second test, and maybe the strongest proof of the HMC method and of our implementation of 
it is a direct comparison with numerical results from a Hamiltonian simulation. The Hamiltonian model 
was defined in terms of spins rather than SU(2) matrix, in order to provide a full proof of equivalence. 
We have carried simulations with both algorithms with the same parameters. We have chosen a 4^^ lattice 
at T 1/20 for several values of the antiferromagnetic coupling to go over the different phases of the 
system. Some of the measures are presented in Table ^ We observe a perfect agreement with precisions 
up to a few per thousand. 

6 Numerical results. 

In this section we present the results of our HMC simulation using the perfect action in the region of the 
Paramagnetic-Ferromagnetic phase transition at vanishing superexchange coupling. We have chosen a 
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fixed temporal length — 40 varying the temperature through a A variation. 

For simplicity on the analysis, we have restricted ourselves to the half-filling case. Due to the hole- 
particle symmetry of the DEM, this can be ensured by setting the chemical potential to zero (see Eg. (p^) ) . 
The study of other band-fillings requires to carefully tune the chemical potential, and will be left for 
further work. 

We have simulated in lattices of spatial sizes L = 4,6, 8, 12, 16 for several values of the temperature. 
We measure every HMC trajectory discarding up to 600 for thermalization (in the worst case). We collect 
between 1000 and 10000 measures at every point. We display our results for the spin magnetization 
(squared) in Fig. |^. The time needed for a trajectory in a 500 MHz Pentium III is about six minutes for 
a 12^ lattice in the critical region with 25 leap frog steps of size 0.02. 

Let us first define the measured observables, and show their general temperature and lattice-size 
evolution, and then later consider in detail their behavior close to the critical region, and measure the 
critical exponents. 

The observables are best defined in terms of the correlation function (the (•) stands for Boltzmann 
average, 

y 

and its Fourier transform, G{k). Then the susceptibility is proportional to the squared magnetization: 

X = G{k = 0) = V{M^) (69) 

It is also very useful to consider a finite-lattice correlation-length, in terms of the minimum allowed 
momentum in a finite lattice ||l^, fcmin = (27r/L,0,0) 

£ = ^- ( l\ (70) 

2sin(fc^in) \G{k^,n) J 

Notice that the above definitions use non-connected correlation functions, and therefore the above 
correlation-length diverges in the ferromagnetic phase like 0{L^^^/^). In the thermodynamic limit, 
£ diverges at the critical point like It]"", (t is the reduced temperature). The critical behavior for x is 

In Fig. H we show the temperature and lattice size evolution of JVP. There are several features to be 
noted. The first one is that the behavior of the L = 4 lattice is rather pathological. We believe that this 
evidences better than any other example the need for larger lattices simulations of spin-fermion models. 
It is also interesting to notice the larger lattices rapidly tend to their thermodynamical limit, out from 
the critical region. Finally, we observe that the low temperature behavior of is linear. This can be 
readily understood if we set that the average direction of the magnetization is, say, the third axis. In 
that case, 

M'^1-IpJ2 (('^^)' + (<^'»)') + '^(('^')'' (<^')'' (^V^)') • (71) 

X 

Since the deviations from the perfect ferromagnetic order are proportional to the mean value of a quadratic 
operator, the linear behavior with temperature follows from the equipartition principle, that holds for 
our classical spins at low temperatures. 

In Fig. ^ we show the correlation-length in units of the lattice size. Notice that the curves for the 
different lattices cross at a temperature growing with growing lattice-size. Eventually the crossings should 
occur at the critical point, as dictated by the Finite-Size Scaling Ansatz (see next section). One can also 
observe that £/L is a growing function of the lattice size in the ferromagnetic phase, as it should be. 



6.1 Critical Exponents 

The main question of interest is whether the DEM presents a second order phase transition between 
the paramagnetic phase and the ferromagnetic one, at finite temperature. If the answer is positive, one 
may also wonder about the Universality Class of this phase transition ||2^. In principle, one of the two 
following scenarios should hold: 
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Figure 2: Magnetization squared as a function of the temperature at Jaf —0,p — 0.5 for several lattice 
sizes with Lr — 40. The four (three) leftmost points for L = 6(L = 8) have been obtained with A = 0.25 
and Lr = 80, 160, 320, 640 [Lr = 80, 160, 320). The continuum line is a fit to A(Tc - Tf^, with /3 = 0.37, 
and Tc taken from the L = 8, 12 lattices pair (see text). 



1. The ferromagnetic Double- Exchange interaction is long-ranged enough to enforce Mean-Field be- 
havior |^l[]. The critical exponents would he v = 0.5 and rj = Q. 

2. The interaction is not long-ranged enough: the physical behavior should be the one of the classical 
Heisenberg model in three dimensions The critical exponents would be = 0.71(1) and 
77 = 0.041(2) ||. 

In order to decide which of the above possibilities hold, we have applied the quotients-method 
to the Finite-Size Scaling Ansatz |Q. We recall briefly the basis of this method. Let O be a quantity 
diverging in the thermodynamical limit as t~^° {t = T /Tc — 1 being the reduced temperature). We can 
write the dependence of O on L and t in the following way [M 



0{L,t) = L^"/" 



L 



(72) 



where Fq is a (smooth) scaling function and (— oj) is the corrections-to-scaling exponent (e.g., —lu is the 
leading irrelevant exponent of the Renormalization Group transformation). This expression contains the 
not directly measurable term ^(oo, t), but if we have a good definition of the correlation length in a finite 
box ^(L, t), Eq. (p2h can be transformed into 



0{L,t) = 1"=°^" 



Gc 



L 



(73) 



where Go is a smooth function related with Fq and F^ and the term has been neglected because we 
are simulating deep in the scaling region. We consider the quotient of measures taken in lattices L and 
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Figure 3: Correlation-length in units of the lattice size as a function of the temperature at Jaf — 0, 
p = 0.5 for several lattice sizes. 

sL at the same temperature 

Then, the main formula of the quotient method is 

Qo|Q,., = s-«/'^ + 0(i^-") , (75) 

i.e., we compute the reduced temperature t, at which the correlation length verifies ^{sL,t)/^{L,t) = s 
and then the quotient between 0{sL, t) and 0{L, t). In particular, we apply formula ( [75| ) to the overlap 
susceptibility, x, and the /3-derivative of the correlation length 9tCi whose associated exponents are: 

Xdri = 1 + 7 (76) 
- {2-r,)v. (77) 

Notice that Qo\q^=s '^^^ be measured with great accuracy because of the large statistical correlation 
between Qo and Q^. It is also very important that in order to use Eq. ( |75| ) one does not need the 
infinite- volume extrapolation for the critical temperature. 

In practice, what we do, is to perform a cubic polynomial fit to ^/L as a function of T, on the critical 
region and use the obtained continuous function on the quotients formula (|75|). We find 

1^6,12 = 0.75(4), = G.1284(9)t (78) 

1/8,12 - 0.72(9), Tc = 0.1379(6)t . (79) 

The above results are certainly compatible with the classical Heisenberg model exponent, 0.71(1), and are 
2.5 standard deviations away from the Mean-Field result, 0.5. The estimate of the critical temperature, 
shows a considerable lattice size dependency (it can be shown that the crossing point tends to the critical 
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Figure 4: Finite Scaling behavior of xl^'^^'^ ^ a-s a function of ^/L, notice the data collapse. The arrow 
signals the value of ^/L at the critical point. 

point as J^-^I^-^ ^ uj being the universal scaling-corrections critical exponent [^). Using the crossing 
point for (8,12) as an estimation of the critical temperature, we can perform a fit of the magnetization 
squared to the function A{T^ - Tfl^. In Fig. | we show a fit with the 0(3) exponent (3 = 0.37||2) (solid 
line). The MF value would correspond to a linear behavior {(3 — 0.5). It seems therefore safe to conclude 
that the second scenario is the one realized in Double-Exchange materials with continous transitions, 
which should have non MF critical beahaviour. Let us however remark that a really accurate measure of 
critical exponent would require the extension of the reweighting techniques p5|] to these models. 

It is amusing to observe that the ratio between the real critical temperature at half filling, Tc w 0.14f, 
and the variational Mean Field estimate, = 0.19t 0, is rather similar to the corresponding ratio for 
the three dimensional classical Heisenberg model (Tc — 1.443Jaf = 2Jaf)- 

We finally perform the plot suggested by Eq. (fz^): x/L''/'^ should be an universal function of i/L. 
This seems to be rather well satisfied by our data, with the critical exponents 7 and v of the classical 
Heisenberg model in three-dimensions. 



7 Conclusions and Outlook 

We have proposed a general numerical method for studying systems consisting of classical degrees of 
freedom coupled to fermionic fields. The method is based in the Path Integral formulation of Quantum 
Mechanics that allows to work in a classical space-time lattice where powerful Monte Carlo techniques, 
as the Hybrid Monte Carlo method, are applicable since no sign problem arises. 

As an example, we have describe explicitly the formulation of the method in the case of the Double 
Exchange Model, observing that is convenient to use a mapping of the spin to SU(2) matrices to avoid 
singularities related with the parameterization of the Berry phase. 

We have also shown that when the spectrum of the single-particle Hamiltonian matrix is bounded, 
it is possible to work directly in the continuum-time limit, using a perfect action thus avoiding the need 
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for a Trotter-Suzuki extrapolation. We have also shown how to eliminate the spurious dynamical effects 
induced by the empty fermion system, when the spectrum is unbounded. 

We have finally presented some numerical results. First we have described some consistency checks 
and then we have studied a property of the model with direct physical interest as the Paramagnetic- 
Ferromagnetic transition. We have studied the phase transition at half-filling, where the transition 
temperature is highest, for simplicity. We have shown that the Finite Size Scaling Ansatz is well satisfied 
for this model. The critical exponents have turned out to be fully compatible with the ones of the three 
dimensional classical Heisenberg model, and incompatible with the Mean-Field prediction, as expected on 
Universality grounds if the interactions are not extremely long-ranged. This conclusion was definitively 
out of reach with the lattices that could be simulated with previous methods. 

Work is in progress for the study of the full phase diagram of the model, (p, T, Jaf)- We are also 
planning to use this Monte Carlo method for the study of models with several electron orbitals and/or 
phonons. 
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A Proof of Eqs. and 

We recall that l/T = L^-A and that the sums (or products) in po run over the Matsubara frequencies 
We apply the Poisson summation formula llq| (valid for a continuous (27r)-periodic function) 

^ E /(^^o) = (-1)^ f (80) 



to the RHS of Eq. (M§): 



s — -\-oo 



±y 2r = y (-i)W^e--* (81) 

Po S — — 00 ^ 

s= — oo \^\ — ^ 

where the orientation of the contour is positive. For s < it is useful to perform the integration in 
w = 1/z. 

When n > E, only the terms s > contribute, and one obtains 

J_ V — - V [-e^^^(^-^)l ^ = - f83l 

Lr ^ ei» - e^(s-p) \ 1 , ' 

PO 5=0 i + C -1 

while if /i < i?, we need to consider only s < ~1 arriving to 
1 e'Po 



V — = - V [-e-^^-f^-^'l ' = - (84) 



po s=l 



To prove the relation (^^ , we start noting that (for L^. even) the products in its LHS can be grouped 
in pairs of nonzero complex conjugates, so it is possible to write 

■Q (c^-^+'Po _ e^^^) = e«(''^^'^) (85) 

Po 



15 



where the function G(/i, A, E) is real. To obtain G we first compute the fj. derivative 

dG _ ^ Ae^^+'P° 

= ^^-i + e^i(^-A^) ^^^^ 
= Alog(l + e-^^^(^-^)). (88) 

From this relation, we know G{fi, A, E) up to a /i-independent term Go(A, £') 

G(/^, A, E) = log (l + e-^^-(^-^)) + Go(A, E) . (89) 

To evaluate Go it is enough to observe that 

lim e'='('^'^^-^) = (-l)^-e^-^--^ (90) 

fl—^ — QO 

lim logfl + e-^^-^^-f^A = 0, (91) 
consequently, Gq ~ E/T and Eq. ( |39| ) follows. 

B Integrals over SU(2) and the sphere 

In this appendix, we want to show that a generic integral over the sphere 



S2 



27r 



D<j>fi$) = j-I dip I sinOde fie,ip), (92) 



can be substituted by an integral over the SU(2) group (with Haar's invariant measure). 

In order to see how can this be possible, we start noticing that, without loose of generality, the function 
depending on the vector variable, /(</>), can be considered as a function of the matrix (0 • <?), because 



.iTr 
2 



, i = 1,2,3. (93) 

Now, one can always find an SU(2) matrix [/[cj)], such that 

Um$-<j)U^$]^a3. (94) 
An explicit choice is given in Eq. (^. There are two important facts to be noticed: 

Two SU(2) matrices, V and W verify V^^asV = W^'asW if, and only if, V = e'"'"^!^ for some a, 



-TT < a < TT. 



• For any SU(2) matrix, W, there is a point on the sphere 4>w, such that W^'cr^W — (pw ■ a 
Therefore, the SU(2) group can be parametrized as 

W = e'"'"-^C/[(^] , -TT < a < TT , = 1 . (95) 
The above considerations lead us to the following chain of equalities: 

Ai>!($) = / ^4>!{4>-a) (96) 

S2 JS2 



d(?/([/t[0]a3C/M) (97) 
'^^h f ^ (c/^[0]e^'"^"a3e-^"C/[(?]) (98) 
dW fiW^'asW) . (99) 



5(7(2) 
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So we see that there is at least one integration measure over the SU(2) group, for which our objective 
can be accompUshed. The only thing that still remains to be done is to show that the above integration 
measure is the proper Haar measure. It will be convenient to recall that the Haar measure is the only 
one which is right invariant [pl]| , namely for any function F over SU(2), and any SU(2) element V, one 
should have 

/ dWF{W)^ I dWF{WV). (100) 

JSU{2) JSU{2) 

But, it is easy to see that iiW = e''^3"C/[(^], then 

U[$]V ^e^^'^'^^''^U[Rv$], (101) 

where Ry is the S0(3) rotation matrix associated with the SU(2) matrix V, in the canonical homomor- 
phism between both groups 

[Rv$]-& = V^{$-a)V . (102) 

At this point wc can just go downhill: 

/ AWF{WV) = / d,^-!- / rfa Ffe'('^(^''?)+")'^^C/[i?y. 
Jsu(2) Js-^ 27r V 



(103) 

(e'""^[/[i?y0]) (104) 

/ f da F (e'^^'^Ui^]) . (105) 

Js2 2n \ I 



>^ I da Fie" 



In the above expressions, the second equality follows from the periodicity in a of the integrand, while the 
third is a consequence of the rotational invariance of the measure on the sphere. 

In order to formulate the Molecular Dynamics equations of motion, one needs to know how to calculate 
derivatives on the SU(2) group. For the shake of completeness, we give here the pertinent definitions, 
but refer to |l^ for a complete exposition. 

One defines three different derivatives over SU(2) (one per group generator) 



(106) 



If / is a smooth function of the matrix element of [/, Ua.p, we have 



UUa.B 



If it depends in the full lattice configuration, {Ux\, we define 



5.,/(f/) = E Jtt^ • (108) 



C The exponential of the single-particle matrix 

In this appendix, we show how to numerically deal with the exponential of a matrix, like the Double- 
Exchange single-particle Hamiltonian matrix, with eigenvalues verifying —6 < En < 6. 

Let us call the coefficients of the Legendre polynomials expansion of the fmiction e^^^ for a: G [—1,1]. 
We can write 

oc 

gAi/oEM = C^P„(iJDEM/6) . (109) 

n=0 
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In the following, we shall use the shortcut H = -ffDEM/6- In practice we use the truncation 



N 



n=0 



that correspond to a Hamiltonian 



log Q{N,X) 
X 



The truncation error is quantified through the function 
Rn[x,X) = 



A 



X e [—6, 6] 



(110) 



(111) 



(112) 



that would be zero if the real exponential was calculated. For instance, Rio{x, 1/2) < 2 x 10^*^ for all the 
interval. 

To preserve the numerical stability is better to use the recurrence-relations of the Legendre polynomials 
than their actual expressions in terms of H. Starting from 



we will use (for ti > 1) 



PaiH) = l, Pi{H)^H, 



Pn+i{H)\v) = ^-HP,,{H)\v) P,^_,{H)\v) 



n+1 



n + 1 



(113) 



(114) 



Notice that since matrix H is sparse (6 non- vanishing matrix element per row) , the truncated expression 
for the exponential can be calculated in order V operations. 

In the HMC, to integrate the equations of motion, we need to know the matrix elements 



N 



.SPn{H) 



(115) 



From (114) we can write a recursive relation for the derivative. However it would mean a recursion 
(involving 0{V) multiplications) for each lattice site. This would make a total of 0{V^) operations. 

Fortunately, it is possible to obtain the matrix elements with 0{V) operations. To this end we use 
the double expansion 



SPnjH) 



71— In — 1 — mi rrr 

= E E L(-l^.PmAH)^PrnAH). 
mi— 1712—0 ^ 



(116) 



In this equation Lm\.m2 £^re symmetric in mi,TO2 and vanish for mi + m2 > n. They can be obtained 
from the following relations 



• If mi + m2 < n — 2 



(„+i) _ 2n + l 



2mi - i""'"!-!^'"^ "^27711+3 "1+1^™^ 



mi („) 



n + 1 



^(n-l) 

7711 .TTlo ' 



• if n — 1 < 7711 + 7712 < 77, with 7771 7^ 

(77+1) _ 277 + 1 7771 



Li 



-L 



(77) 



• Finally 



77 + 1 277;,l — 1 
(77+1) _ 277 +1 



7771 — 1,7772 ' 



77+1 



In terms of the L coefficients we can write 



(Glflc'J-^ \F) = ( E (GlPrrniH)) §- (e " E'" PrnAHm] 



(117) 

(118) 
(119) 

(120) 



77 = 



^ 7771=0 



^77 = 



18 



References 

[1] E. D. WoUan and W. C. Koehler, Phys. Rev. 100, 545 (1955); D. I. Khomskii and G. Sawatzky, 
Solid State Commun. 102, 87 (1997); J. M. D. Coe, M. Viret and S. von Molnar, Adv. in Phys. 48, 
167 (1999). 

[2] C. Zener, Phys. Rev. 82, 403 (1951); P. W. Anderson and H. Hasegawa, Phys. Rev. 100, 675 (1955). 
[3] A.J. Millis, B.I. Shraiman and R. MuUer, Phys. Rev. Lett. 77, 175 (1996). 

[4] M. J. Calderon, J. A. Verges and L. Brey, Phys. Rev. B59, 4170 (1999), and references therein. 

[5] A. Moreo, M. Mayr, A. Feiguin, S. Yunoki and E. Dagotto, Phys. Rev. Lett., 84, 5568 (2000), and 
references therein. 

[6] R.T. Scalettar, D.J. Scalapino y R.L. Sugar, Phys. Rev. B34, 7911 (1986); S. Duane, A.D. Kennedy, 
B.J. Pendleton y D. Rowcth, Phys. Lett. B195, 216 (1987). 

[7] J. L. Alonso, L. A. Fernandez, F. Guinea, V. Lahena and V. Martfn-Mayor, |cond-mat /000347^ , and 
references therein. 

[8] M. Uehara, S. Mori, C. H. Chen and S.-W. Cheong, Nature 399, 560 (1999); J. M. De Teresa, C. 
Ritter, M. R. Ibarra, P. A. Algarabel, J. L. Garci'a-Muhoz, J. Blasco, J. Garcia and C. Marquina, 
Phys. Rev. B56, 3317 (1997); J. Fontcuberta, B. Martinez, A. Seffar, S. Piiiol, J. L. Garcia-Munoz 
and X. Obradors, Phys. Rev. Lett. 76, 1122 (1996); P. Dai, H. Y. Hwang, C. Kloc and S.-W. Cheong, 
Phys. Rev. Lett. 80, 4012 (1998); J. Mira, J. Rivas, F. RivaduUa, C. Vazquez- Vazquez, and M.A. 
Lopez-Quintela, Phys. Rev. B 60, 2998 (1999). 

[9] J.L. Alonso, Ph. Boucaud, V. Martfn-Mayor and A. J. van der Sijs; Phys. Rev. D61, 034501 (2000). 

[10] See e.g. I. Montvay and G.Miinster, Quantum Fields on a Lattice, Cambridge University press, 
Cambridge (1994) 

[11] See for instance M. Creutz, Quark and Gluons on a Lattice, Cambridge University press, Cambridge 
(1983) 

[12] See e.g. T. Inui, Y. Tanabe and Y. Onodera, Group Theory and Rs aplications in Physics, Springer- 
Verlag, Berlin Heidelberg (1990). 

[13] See e.g. M. Imada, A. Fujimori and Y. Tokura Rev. Mod. Phys. 70, 1040 (1998) 

[14] P. Hasenfratz, F. Karsch, Phys. Lett.. B125, 308 (1983); J.B. Kogut, H. Matsuoka, M. Stone, H.W. 
Wyld, S. Shenker, J. Shigemitsum D.K. Sinclair, Nucl Phys. B225[FS9], 93 (1983). 

[15] G. Batrouni et al, Phys. Rev. D32, 2736 (1986); S. Gottlieb, W. Liu, D. Toussaint y R.L. Sugar, 
Phys. Rev. D35, 2531 (1987). 

[16] R. Gupta et al, Phys. Rev. D40, 2072 (1989). 

[17] See, e.g., chapter seven in [|lO|. 

[18] See, e.g., chapter four in JlOf . 

[19] F. Cooper, B. Freedman and D. Preston, Nucl Phys. B210, 210 (1989). 

[20] See e.g. G. Parisi, Statistical Field Theory, Addison Wesley, New York, (1988). 

[21] For recent work on the crossover between non trivial and trivial critical behavior, see E. Luijten, 
Phys. Rev. E59 (1999) 4997, and references therein. 

[22] H. G. Ballesteros, L. A. Fernandez, V. Martin-Mayor and A. Muiioz Sudupe, Phys. Lett. B387, 125 
(1996) 



19 



[23] H. G. Ballesteros, L. A. Fernandez, V. Martin-Mayor and A. Muiioz Sudupe, Phys. Lett. B378, 207 
(1996); Nucl. Phys. B483, 707 (1997). 

[24] See for instance M.N. Barber in Phase Transitions and Critical Phenomena, vol. 8, edited by C. 
Domb and J.L. Lebowitz (Academic Press, London, 1983). 

[25] M. Falcioni, E. Marinari, M. L. Paciello, G. Parisi and B. Taglienti, Phys. Lett. BIOS, 331 (1982); 
A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 61, 2635 (1988). 



20 



